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Abstract 

The Lattice Field Theory approach to the statistical mechanics of a classical 
Coulomb gas [R.D. Coalson and A. Duncan, J. Chem. Phys. 97, 5653 (1992)] is 
generalized to include charged polymer chains. Saddle-point analysis is done on the 
functional integral representing the partition function of the full system. Mean-field 
level analysis requires extremization of a real-valued functional which possesses a 
single minimum, thus guaranteeing a unique solution. The full mean-field equations 
for such a coupled system are derived, as well as the leading (one-loop) fluctuation 
corrections. Two different numerical real-space lattice procedures are developed 
to implement the generalized theory; these are applied to the problem of a charged 
polymer confined to a spherical cavity in an electrolyte solution. The results provide 
new insight into the physics of confined polyelectrolytes. 



1 Introduction 

In a recent series of papers, the application of lattice field theory (LFT) techniques 
[H, H, 1^, §, 1^ to systems of mobile polar ions interacting with immobile macroions 
has led to a calculational framework for extracting free energies, ionic concentra- 
tions, electrostatic potentials, and other physical properties. This formalism, which 
relies on a Hubbard-Stratonovich transformation of a functional integral represen- 
tation of the grand canonical partition function, has a number of antecedents in the 
literature [|, 0, §. In the generalization introduced in |jl|, we were able to obtain 
accurate results at the mean-field level, with precisely computable corrections indi- 
cating the level of accuracy of the mean-field results, for systems of highly variable 
ionic concentration and macromolecular charge. Paper ^ focussed on interactions 
between charged spherical colloid particles (polyballs) in electrolyte solutions. In a 
subsequent paper, a liquid phase assembly of many interacting polyballs was studied 
[0] . The formalism of [H] has been extended to include effects of a spatially variable 
dielectric profile |^ , the finite size of the simple mobile ions forming the electrolyte 
[|], and the presence of mobile dipoles in the gas (solution) p. The objective of the 
present paper is to extend the LFT approach to compute the statistical mechanics 
of a gas of mobile charged ions to a system in which mobile charged polymer chains 
are also present and interacting electrostatically with the mobile ions. 

In Section 2 we briefly review the functional integral representation (and lattice 
field theory discretization thereof) for the grand canonical partition function intro- 
duced in |jl[|. In Section 3 the formalism is extended to include charged polymer 
chains. The mean-field level solution of the full system (ions plus polymers) is then 
presented in Section 4. Our results are similar to, but do not fully agree with some 
recent work of Orland et al [^]. In Section 5 we justify the contour deformation 
introduced in Section 4 to extract the leading saddle-point (or mean-field theory) 
equations of the full ion-polymer system by deriving a general convexity theorem 
for the mean-field free energy. An explicit expression is also derived for the leading 



post-mean-field corrections due to the fluctuations around the saddle-point flelds. 
This means that a precise estimate of the error in the mean-fleld theory is in principle 
possible even in systems of coupled ions and charged polymer chains. Some simple 
illustrative applications of the mean-fleld results of Section 4 to systems possessing 
spherical symmetry, for which the problem can be reduced to the minimization of a 
one-dimensional functional, are described in Section 6. The same systems are treated 
in Section 7 by solving the mean-field equations on a three-dimensional lattice, an 
approach which can be used to study systems of arbitrary shape and symmetry. 
Section 8 summarizes our results and indicates areas for further research. 



2 Functional formalism and lattice field theory for ionic sys- 
tems 

A simple functional integral representation for the grand canonical partition function 
of a system of mobile ions with electric charge density p{r) is obtained by rewriting 
the electrostatic potential energy of the system via a standard functional identity: 

^-^ I 'irdr"^0^ = C f Dxe^-f''^'"''^'-f^^'^'^'^''' (1) 

Here C is an irrelevant constant (independent of p) which will be neglected from 



now on, and the charges are immersed in a medium of dielectric constant e [|10| at 
inverse temperature /3 = 1/kT. For simplicity, assume the ions to be of two varieties 
only: n^ ions of charge +e at locations r ^ and n_ ions of charge — e at locations 
fj' (with 1 < k < n+,1 < I < n^). Of course, if there are no other charged entities 
present, n+ = n^, but for the time being (anticipating the inclusion of charged 
polymer chains below) we leave n± free. The charge density is then 

p(f) =e^(5(f-f+)-e^(5(f-f^"") (2) 

k I 

The effect of the transformation induced by (1) is twofold. First, the nonlocal 
Coulomb potential (which makes the direct simulation of such systems by molecu- 
lar dynamics extremely difficult) has been replaced by a local Laplacian operator. 
Secondly, the nonlinear pairwise coupling of the charged ions has been replaced by 
a linear coupling term where the ions interact only with the auxiliary field. Conse- 
quently, for a fixed x field, the statistical mechanics of the ionic gas is elementary. 
Introducing chemical potentials fi± for the positive and negative ions, the sum of 
the factor exp [i J x{'^p{'^dr) over ion numbers and positions can be explicitly per- 
formed: 






= expie^^^ J^e^'^ + e^'^- J^e-''^) (3) 

where X± are thermal deBroghe wavelengths for the ions. Inserting this result in the 
complete expression for the Boltzmann weight (1), the full grand canonical partition 
function becomes 

where we have rescaled the auxiliary field x ~^ f^X- The generalization of this expres- 
sion to include stationary charged macroions |1|, external single particle potentials 
acting on the mobile ions, spatially varying dielectric constant 0, finite size effects 
[^], and effects due to mobile multipolar entities ^ have been discussed elsewhere 
and will not be treated here. For example, the inclusion of a single particle potential 
V{r) acting on the ions can be effected simply by the replacement e^*^^^ — >■ 6=*=*^^^"^ 
in (^). Such a potential can be used (for example) to exclude the ions from certain 
regions of space. 

The functional integration indicated in @ may be made well-defined by discretiz- 
ing the system on a cubical lattice defined by lattice points r = ai{nx, Uy, n^), where 
ai is a lattice spacing taken small in comparison to the length scales over which the 
field X varies substantially, and n^^ny^rtz are integers, < nx,ny,nz < L — 1, for 
a LxLxL lattice. We shall restrict our calculations to systems with zero total elec- 
tric charge, where the fixed charges are screened by the mobile ones and the fields 
fall exponentially far from the sources, so finite volume effects go to zero rapidly 
as the lattice size Lai is increased. With these notations, the lattice version of the 
transformed partition sum (^ becomes 

^gc = y n ^Xn exp < - X] Xn^nniXm 
n \ nrn 

+ 7+^e*^^>^" +7_^e-^"'^^4 (5) 

n n J 

Here A^^ is the discrete lattice Laplacian operator, which may be regarded as a 
Lr'xL^ matrix (see 0). The dimensionless constants a = aiPe/Arc, '~f± = e^'^^{ai/X±)^ 



have been introduced to simplify the notation. In (5), A now represents the dis- 
crete lattice gradient (nearest-neighbor difference operator), with e the electronic 
charge and ai the lattice spacing. In Section 4 we shall return to the evaluation of 
the integral (5) by saddle-point techniques, after the additional terms implementing 
charged polymer chains are included. 



3 Functional formalism for charged polymer chains 



The path integral representation of a polymer chain ||Tl|, is well-known. We shall 
show below that for a system consisting of many monomers the leading behavior 
is determined simply by the total length M (= total number of monomers) of all 
chains and not by the connectivity or number of the separate chains. Thus the 
configuration of the polymer chain(s) can be represented by a single function x(s), 
where the dimensionless path length variable < s < M . If the typical monomer 
separation is a^ (the Kuhn length Jill), ^^^ partition sum for such a system may be 
written 



0p(r) = / ds5{f — x{s)) (7) 

The field (/)p is the spatial monomer density and the term involving A / ^pdr* prevents 
the monomers from overlapping with each other, A > being a measure of monomer 
excluded volume. As discussed at the end of Section 4, A ~ -B2 — cr^, where B2 is 
the second virial coefficient, which goes approximately like the monomer volume 
(i.e. 0" is the effective hard-sphere radius of a monomer) [l^. For simplicity we 
shall assume initially that the polymer charge density is uniform along the polymer 
chains, with each monomer carrying a net charge of —pe (the generalization to 
varying monomer charge would simply result in an effective Schrodinger Hamiltonian 
with a time-dependent potential (cf. discussion below)). The negative sign indicates 
that we have taken the polymer chain to be negatively charged, as a consequence 
of dissociation of positively charged counterions. This means that (2) for the total 
charge density now becomes [0 



Pir) =e^(5(f-f^)-e^(5( 



r — r I j — pe(pp[r) [6) 

k I 

Repeating the argument leading from (2) to (4) we now find the following ex- 



pression for the full grand canonical partition function of the system 



Zgc = J DxD(f)pDx{s)e' 

The integration over polymer configurations can now be replaced by an equivalent 
Schrodinger Hamiltonian problem. To facilitate this, linearize the dependence on 
the monomer density 0p by introducing a second auxiliary field 

g-t/<^^rfr^ f JJ^^-^Jc.^dr-^xJu;(r)M^df ^^q^ 

Introducing this representation into (9) and using (7), one finds 
with c± = e'^'^^/A^, and 

^ / dsx (s)—ipe3 j dsY(x(s))—iX j dsLo(x(s)) 

defines a functional which is just the Feynman path integral for the imaginary time 
evolution of a particle of mass S/a^ in an imaginary potential i{pe(3x + Acj). In 
fact, we shall show below that the functional integral (11) over x ^"^^ ^ can be 
rerouted through a complex saddle-point at x = —iXc and uj = —iUc where Xc 
and Uc are purely real, so that the evaluation of Zschrix^^) at the saddle reduces 
to a completely conventional three-dimensional Schrodinger Hamiltonian problem — 
namely, the computation of matrix elements of e~^"^ where the Euclidean time 
extent of the evolution is just T = M and 

H = --f V^ + Xco.ir) + (3pexc{^ (13) 

o 

For large M, any matrix element of e~^'^ has the asymptotic form 

log ((... y^'\ ...)) ^ -ME^ + 0(1) 
where Eq is the ground state eigenvalue of the Hamiltonian H |15|. Different bound- 



ary conditions (periodic, open, etc) imposed at the ends of the polymer chain (or 



even a finite but fixed number of splits in the chain) correspond to the choice of 
different initial and final states in the above formula, contribute to the 0(1) edge 
correction, and are subdominant for large M . 
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4 Mean-field (saddle-point approximation) theory of cou- 
pled ionic-polymer system 

The equations determining the saddle-point configuration fields Xo ^c are obtained 
by setting the variational derivative of the exponent in the full functional integral 
(11) to zero. After rotating the fields to the negative imaginary axis (the need for 
this rotation will be justified below when we discuss the fiuctuation corrections), 
this exponent becomes 

F = Jdfl^ jVxcf + ^uj! + c+e^'""' + c.e-P'A - MEo{Xc,cOc) (14) 
The functional derivatives determining the saddle-point solution are then 

= /5pe|^o(r)P (15) 



^^0 .__... ._-,M2 



^Xc(r) 

A|^o(rir (16) 



SEo .,.,,^^,2 



6uJc{r) 
^ ^^ ^ -V\c{r) + c+e^^^^("^ - c.e-^'^^^"^ - Mp\^o{r)\^ = (17) 



Pe6xc{^ 47re 

1 SF 

Here the wavefunction \E'o(^) is assumed unit-normalized: 

'"|^o(^)Pt/f = 1 (19) 

Using (18), the auxiliary field Uc may be eliminated completely, leaving the pair of 
coupled nonlinear equations as the complete mean-field solution to the full interact- 
ing ion-charged polymer problem: 

^^\{f^ = c+e^'^^=('')-c_e-^"^^('^-Mj9|^o(r)|' (20) 

47re 

-fV2^o(r) = AM|^or^o(r)+/9j9eXc(r)^o(r)-^o^o(r) (21) 

o 

As mentioned previously, the inclusion of single particle potentials, which can be 
used to enforce exclusion regions for either the ions or the monomers, is straightfor- 
ward. Using the notation Vc{r) (resp. Vm{r)) for potential fields acting on the ions 
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(resp. monomers) (20-21) become 

-V\cir} = c+e'^^^=("^"^=^'^'^-c_e-''^^=^'^'^-^=^'^'^-Mp|^o(r)|' (22) 



-^V'^o(r) = AM|^or^o(r)+/?peXc(r)^o(r)-(^o-V;«(r))^o(r) (23) 

Recalling that the parameters c± are exponentials of the chemical potentials fi± 
for positively and negatively charged ions, the numbers of these ions must be fixed 
by suitably adjusting c±, using the easily derived relations: 

n± = c^^MM = c^f e^^^^^-^'^dr (24) 

oc± J 

while the number density of monomers is given by M|\l'o(^) P which clearly integrates 
to the total number of monomers M. Charge neutrality follows from integrating (22) 
over all space (correct either for periodic boundary conditions or for an electric field 
falling to zero at the system boundaries): 

= /"v^Xc(r)rff ^ = n+ -n_ -Mp (25) 

=^ n+e = n_e + Mpe (26) 

For systems of spherical symmetry, the equations (22-23) reduce to a pair of 
coupled nonlinear ODEs: 

r^'lL = ^ (^^e^M/'-^^M - e-e-^«/^"^^M) - ^^(r)^ (27) 

dr ^ / f 

\^ = ^Fir)Gir) + ^Girr-iEo-VUr))Gir) (28) 

where the rescaled radial functions F, G are defined as 

F{r) = I3erxc{r) (29) 



G{r) = yAMr^o(r) (30) 

The radial variable in (27-28) is measured in units of the monomer size Op, and we 

have defined dimensionless variables 

Ae 

47r/3pe^a^ 

e± = ^ (32) 

P 
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^ ~ A-0^,2n2 ^'^^^ 



In terms of these new variables the charge and number density constraints (24-26) 
become 



,,/,.V^-«.,. . ^n, (33) 



A 

A 



47rpap 
Girfdr = ^M 
A 



47ra3' 



3 ,n^-n_) (34) 

In fact, we shall show below that the solution of the coupled nonlinear equations 
(27-28) is most readily accomplished by returning to the free energy expression 
(14), which represents a functional of two scalar fields Xc, ^c- This functional will 
be shown to be convex and bounded below, with a unique minimum at the field 
values satisfying (27-28). 

For non-spherically-symmetric systems, a practical approach is again provided by 
the lattice field theory discretization of [|l| . The inclusion of the polymer terms in 
(5) is quite straightforward — discretizing (11), we obtain (setting the single particle 
potentials Vc, Vm to zero for simplicity) 

^gc = y n dXndUH exp < - XI Xn^nrnXm " -af JZ ^^ 
n \ nm ft 

+ 7+ E e'^^"^' + 7- E e""'''^'^ - MEo(Xn, ^n) I (35) 

n fi ) 

where Eq is the ground state energy of the associated discrete Hamiltonian matrix 



a" 



Hfnn = -ITl^rhn + i^^n + ^/^Pe^n (36) 



6a? 



Of course, at the complex saddle-point of the integral (35), x = —^Xc and uj = —iujc 
with Xc) ^c real fields, so the Hamiltonian matrix (36) is real symmetric, with well- 
defined real ground state energy. After the rotation to the imaginary axis, we find 
the following expression for the discretized free energy functional (corresponding to 
the continuum expression given in (14)): 

^= -f EXrnA^nXn+^«fE^i + 7+Ee^'^'+7-Ee-^''''^-Mi^0(Xn,C^n) (37) 
fhn n n ft 
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where we have dropped the subscript c indicating the saddle-point (mean-field) value 
for the fields to avoid notational overload. Now the ground state energy, as promised, 
is a perfectly real number, namely the lowest eigenvalue of the matrix 



a? 



Hrrin = --^l^rnn + ^^^n + (^V^^Xn (38) 

6af 
For a long Gaussian polymer chain in a confined region and subjected to an 



external potential Wir)^ it is well-known UTTj] that the ground state of the three 
dimensional Schrodinger operator 

i^ = _5^V2 + /5H^(f) (39) 

6 

determines the equilibrium properties of the polymer chain. In particular, if \E'o(^) 

is the unit-normed ground state eigenfunction associated with Hamiltonian (|39|) 

(with \E'o = on the confining boundary surface) then |\l'o(r)p is the normalized 

monomer density. Equivalently, if there are M monomers in the polymer system 

then M|\E'o(r)P is the number density of monomers in the system. Comparing (|39|) 

to (38), (18) above, we see that in our problem there is an effective potential 

W{f) = fcTAM|^o(r)|^ +pexc(r) 

Since both terms in this potential are functions (explicitly or implicitly) of \E'o) our 
effective potential W is evidently a mean-field potential. Specifically, the electric 
potential Xc depends on the monomer charge density (hence "^o) according to (20). 
The other contribution to the effective potential, i.e. /cTAM|\l'o(r)P, is due to short- 
range excluded volume interactions between monomers. Its origin and form can be 
understood as follows. Let U{r) be a short range repulsive potential via which 
all monomer pairs interact. Then, if there is a distribution nm{r) of monomers 
in the system, the repulsive potential (call this W, suppressing the electrostatic 
contribution) experienced by a test monomer inserted at point f is: 



W{r) = df'U{f-f')n^{f') (40) 

= n^{r) f dfU{f) (41) 
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The second line follows from the first under the assumption that the pair potential 
U is short range compared to the length scale on which the monomer density Um 
varies. Using the connection (noted above) that M|\E'o(r)p = nm{r), we immediately 
identify the parameter A as 

X = (3 f U{f)df (42) 

This can in turn be connected to the second virial coefficient B2, which can be used 



to ascribe an effective "hard sphere" radius a to the monomer ||13|| : 



p / U{f)df 



df 



-PU{r) 



2B2 = 4nay3 



In this way we can connect A with the effective size of a monomer. 



(43) 
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5 Corrections to mean-field tiieory — convexity of the free 
energy 

In the preceding section we claimed that the functional integral (11) is dominated by 
a complex saddle-point where the fields x? ^ take pure imaginary values —iXc —iuJc- 
To justify this assertion, we need to ensure that the deformation of the contour of 
the field integrations is such that the integral passes through a proper saddle-point, 
with a maximum of the real part of the exponent, and that the saddle chosen is 
the dominant one globally. A rigorous discussion presupposes that the functional 
integral (11) has been made well-defined, say by discretization as in (35). As a 
consequence we may be sure that the contour deformation is possible in the first place 
as a result of the analyticity of Zschr = Tre"^'^*^^"''^"^ as a function of the variables 
Xni ^n for an arbitrary finite-dimensional matrix H . We remind the reader that the 
result 2'schr = Tre~^'^ holds for closed polymer chains but that different boundary 
conditions induce subdominant edge effects of order 0(1/M) for long chains with M 
monomers (cf. discussion following (13)). The analyticity of the remaining part of 
the integrand in (35) is obvious. 

For simplicity we shall temporarily return to continuum notation, although the 
reader is warned that ultimately the functional integral being discussed must be 
explicitly defined by a cutoff procedure. Writing 

Xif) = -iXcif) + xif) 

uj{r) = —iujc{r) + oj{f) (44) 

we may expand the integrand of the functional integral (11) keeping terms up to 
second order in the fluctuation flelds x? ^'■ 

Zgc = e / Dx[r)Duj{r)e J V»"' ' ' ^ ^ ^ J 

^ g-M/drdr'(AiD(r)+/3pex{r))Gc(r,r')(A'i{r')+/3pex(r')) U^\ 

where the prefactor e^ contains the entire mean-field result (14) and the "one- 
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loop" fluctuation corrections to mean-field are contained in the remaining Gaussian 
functional integral over the fluctuation fields x, u. The crucial point is that this 
integral is convergent as a consequence of the positivity of the Green's function 
Gc{r, r ') giving the second order variation with %, uj of the ground state eigenvalue 
Eq of the hermitian Hamiltonian (13). Let \I^„(r) be a complete orthonormal set of 
eigenfunctions of H, with the corresponding ordered eigenvalues En, En < Em for 
n < m. Then standard second order perturbation theory gives, as a consequence of 
(44), 






Eoix,^) - Eo{xc,uJc) + J drdr\\uj{f) + l3pex{f))Gc{f,r'){\uj{r') + Ppex{r')) 

^„(f)^„(r") 

through terms of second order in the fluctuation fields. Note that the first order 
perturbation shift in Eq is canceled at the saddle-point by the first order variation 
in the exponent in (11). Furthermore, the positivity (strictly speaking, positive 
semidefiniteness) of Gc is manifest. To see this, note that a necessary and sufficient 
condition for Gc to be positive semidefinite is 



df df'f{r)Gc{f,f')f{f') > (47) 

for any function /(r). Using the definition given in (|46|) , 

[ ,^[ ,^,.,^^ ,^^,..,^,. ^ [/rff^o(r)/(f)^„(f)]' 
dr dr f{r)Gc{r,r )f[r ) = }^ (48) 

■J -I n^Q ^« ~ ^0 

Since the denominator of each term on the r.h.s. is positive and the numerator is 
non-negative the condition for positive semidefiniteness (^7|) is satisfied. 

It should be emphasized that the positivity of the full fiuctuation kernel in (45) 
holds for arbitrary real fields Xa ^c — it is not essential that they satisfy the saddle- 
point equations (20-21) (although we shall of course eventually demand that they 
do!). The value of this observation is that it is equivalent to a statement of con- 
vexity of the free energy functional F{xc-,^c) in (14) for arbitrary values of its field 
arguments, as the exponent in the fiuctuation integral (45) is essentially the second 
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functional derivative of F. The functional F{xc, ^c) is clearly bounded below (as 
the lowest eigenvalue of H in (13) can grow at most linearly with Xc or cUc). If it 
is everywhere convex, it must have a unique minimum at exactly the field values 
Xci^c satisfying (20-21). Thus the stable minimization procedures employed in ^ 
for solving the lattice field theory of systems of polyballs and mobile ions are guar- 
anteed to work here also, in the presence of long charged polymer chains, provided 
an efficient numerical technique is employed! 

We now return to the task of evaluating the Gaussian fluctuation integral (45). 
It is convenient to change the fleld integration variables by replacing the excluded 
volume fleld Cj{r) with the linear combination 

cr(r) = XCj{r) + l3pex{r) (49) 

so that the fluctuation integral in (45) becomes 

^ ^~M J drdr'a{f)Gcir,r')a{r') /ggN 

whence 

Fi = — InDetiT (51) 

where K is the kernel 



K 



-fA + i£|P + i¥(..e''«. + c_e- 



2A 



V -ff MG. + k 



This determinant must be rendered well-deflned by an explicit cutoff procedure, 
such as the lattice. On a lattice with N points, we then have the problem of evalu- 
ating a 2Nx2N determinant. If the polymer is restricted to a subregion of Np points, 
the lower right hand NxN block in K is non-sparse only in a NpxNp subblock of the 
kernel Gc- The evaluation of Gc is facilitated by the observation that 

E '^''^J^'^f^ = hm ({H - E,) \ \ (f, r ') (52) 

with the required inversion involving only a sparse matrix {H — Eq)"^ + rj'^. 
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6 Mean-field results for spherically symmetric systems 

For spherically symmetric systems, the extremization of the functional F in (14) is 
much simplified. It is convenient to express all distances in units of the effective 
monomer size (Kuhn length) a^, and to define rescaled radial functions /(r), g{r) 
as follows 

f{r) = (3exc{r) (53) 

g{r) = r'^Qij-) (54) 

It is also convenient to rescale the auxiliary field ooc to a dimensionless one, h{r) = 
Ucttp. Introducing dimensionless parameters r], ^±, and ( 

e± = 47rc±4 (56) 

C = 47r4 (57) 

Op 

the saddle-point functional can be written as a one-dimensional integral 

F = j\v[^] +^Ch{rf + ^+ef + ^-e-nr^dr-MEo (58) 

where the radial wavefunction g{r) of the associated Schrodinger Hamiltonian (13) 
satisfies 

l^^2 = [-ih{r)+pf{r)-Eo^g{r) (59) 

The rescaled activity coefficients ^± must be constrained by the appropriately rescaled 
versions of (33): 

We have devised the following efficient procedure for the minimization of F in 
(58). The arguments of the preceding section establish the convexity of F, implying 
a unique minimum. After discretizing the r variable, the rescaled electrostatic po- 
tential function /(r) and polymer density function h{r) become finite arrays which 
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can be updated alternately by the following procedure, which at each step takes us 
closer to the unique global extremum of F: 

1. Choose a reasonable starting value for the fields /, h. 

2. Define a new field a{r) = ^h{r) + pf{r), so that Eq in (58) is the lowest 
eigenvalue of the operator H = — i^ + cr(r). Once r is discretized, H becomes 
a tridiagonal matrix. For the rest of the calculation, we minimize with respect 
to /(r) and <j{r). 

3. Minimize F with respect to cr(r) for each discrete value of r. This requires a 
rapid calculation of Eq as a function of clr). The ground state eigenvalue of 
a tridiagonal matrix (indeed, any ordinally located eigenvalue) can be readily 
extracted by Sturm sequence methods [Q, and the functional F minimized 
quickly with respect to a{r) by golden section bracketing |l^. The latter 
approach is foolproof as we have a strictly convex dependence on a{r) for all r. 

4. Minimize F with respect to /(r) for each discrete value of r. This is trivial as 
the dependence of F on /(r) is explicit. 

5. Iterate until F stabilizes at a minimum to some preassigned tolerance and/or 
the saddle-point equations are satisfied to a desired degree of accuracy. 

This algorithm was applied to a system consisting of a charged polymer with 1000 
monomer subunits trapped in a spherical cavity of radius lOOp, with each monomer 
carrying a charge -O.le. The ions (r7,+=200 positive and n_=100 negative ions) are 
free to move in a larger spherical region of radius lOOa^. The parameter rj was 
taken to be unity. The results shown correspond to 1000 iterations, after which 
the free energy is stabilized to 5 significant figures (such a run takes a few minutes 
on a 400 MHz Pentium processor). The effect of variation of the excluded volume 
parameter ( on the monomer density "^oir) is shown in Fig. (1). As for the case 
of the uncharged polymer, increasing the excluded volume parameter results in a 
flattening of the distribution near the origin. However, in the presence of charge, 
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Figure 1: Monomer distribution for varying C, 
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Figure 2: Electric Potential /(r) for varying ^ 
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we now find that (at least for tlie parameter range studied fiere) tlie electrostatic 
repulsion sufficiently counteracts the tendency of the monomers to crowd into the 
central region to produce a dip in the distribution for small r. The mean-field 
rescaled electric potential /(r) is essentially constant over this variation of (, as 
shown in Fig. (2) ||T^. The exponential screening outside the cavity of radius lOop 
confining the polymer is evident. 

The effect of varying monomer charge p (with the excluded volume parameter ( 
held fixed at 20.0) is illustrated in Figs. (3,4). In particular, as p ^ we recover the 
distinctive fiat behavior at small r characteristic of uncharged polymers (Fig. (3)). 
In these plots, the negative ion number is held constant at 100, with the number of 
positive ions adjusted to give charge neutrality. 
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Figure 3: Monomer distribution for varying p 
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Figure 4: Electric Potential /(r) for varying p 

Finally, we may also study the effect of varying salt concentrations. Increasing 
the background ion density results in higher screening of the electrostatic potential, 
so the effect on the monomer distribution is similar to that obtained by varying 
the average monomer charge p. In Figs. (5,6) we show the monomer distribution 
and electrostatic potential /(r) for fixed C=10 and p=0.1, varying the number of 
negative mobile ions n_, with tij^ = ri- + 100 for charge neutrality. 
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7 Mean-field results using three-dimensional lattice field the- 
ory 

In the previous section we showed how the theory developed in this paper can be 
apphed to the problem of a charged polymer confined to a spherical cavity and im- 
mersed in an electrolyte solution. We now present the solution of the same problem 
using three-dimensional (3D) lattice field theory, which does not hinge on special 
symmetry properties of the system, and thus illustrate a numerical procedure for 
dealing with systems of arbitrary shape and complexity. Colloidal suspensions in 
polyelectrolyte solutions, which may be useful for for a variety of technological ap- 
plications, such as new optical materials and devices (e.g., narrow-band optical 
rejection filters, pump-probe laser apparati, optical display panels, etc. [^ pT|), 
are systems of this type. 

We will solve the discretized versions of equations (22) and (23) simultaneously 
on a 3D lattice. It is convenient to multiply equation (22) by af, where ai is the 
lattice spacing, and to rescale according to: f{r)—*(3eXc{^),i^Q{r)—^ai ^/Jo{f). Thus, 
all variables and parameters become dimensionless and the two discretized equations 
are: 

a Y. ^nrnfrn = 7+6^'^-^'^ - ^_e-f^-^« - Mp^jl (61) 

ffi 

-^ J2 ^nm^rh = -^^H^n + Pfn^n " ^O^n (62) 

\)Cli -, CI] 

' m ' 

where 

7± = ^T^, (64) 

and the wavefunction is dimensionless and normalized according to the rule 

n 

We solve equations (61) and (62) simultaneously using the following procedure. 
First, equation (62) is solved for /« = and ignoring the non-linear (monomer 
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repulsion) term. The resulting ipfi (wave-function for a free particle in a sphere) 
is given to equation (61), which is solved at each lattice point by the Newton- 
Raphson method. The process is repeated and the coefficients 7-1- are updated 
after each iteration with the current field, until a predetermined desired accuracy 
is achieved. Then the resulting /^ is fed into equation (62) which is solved for a 
new ipa to be given to equation (61). Equation (62) is solved to a predetermined 
desired accuracy by the Lanczos approach, which is appropriate for a large sparse 
matrix like the one representing our Hamiltonian. This method of computation is 
very well-suited for implementation on massively parallel platforms which should 
make it possible to study even very large lattices with this approach. From then 
on, the potential f^ to be given to (62) for the next iteration is updated slowly by 
adding a small fraction of the new fa to the old one, obtained from the previous 
iteration. The same overrelaxation procedure is used for updating ■?/'| in the non- 
linear term of the Schrodinger equation (62). Such a gradual iteration procedure is 
necessary in order to avoid an unstable bifurcation between two unphysical states 
(which is commonly encountered when solving non-linear differential equations), 
and it converges to the simultaneous solution of the two equations. We have already 
shown that the functional in (14) has an unique minimum, the condition for which 
is given by the two equations, (61) and (62). Therefore, once we have converged 
to a solution of these two equations, we are guaranteed to have reached the unique 
mean-field solution of the problem. 

We have applied the procedure described above to a system of a negatively 
charged polymer of 1000 monomer units, confined to a spherical cavity of radius 
20ap, which is immersed in an electrolyte solution confined to a larger sphere of ra- 
dius 40ap. The Kuhn length Op has been chosen to be SA. To illustrate the stability 
and accuracy of the procedure, we compare some of the 3D results with the ones 
obtained through the one-dimensional (ID) calculation of Section 6, which are prac- 
tically exact. In Fig. (7) we show how the 3D results for V'o(r) approach the exact 
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ID result as the number of the lattice points on the side of the 3D cube contain- 
ing the system, L, is increased from 40 to 60 and 80, for the following parameters: 
( = 15, p = 0.1 and n_ = 6, where n_ is the number of the negative ions in the 
system, while the number of the positive ions is adjusted so that electroneutrality is 
preserved, and the relationship between ( and A (in equation (62)) is given in (57). 
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In Fig. (8) we show the effect of variation of ( on the probability distribution 
ipQ{r). For comparison we have included two results obtained by the ID method of 
Section 6. The results are analogous to the ones in Fig. (1). Similarly, the rescaled 
electrostatic potential /(r) varies little with (, Fig. (9). 
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Figure 9: Electric Potential /(r) for varying C, 
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In Figs. (10) and (11) we present results for ipoi'r) and f{r) for varying monomer 
charge p, fixing (^ = 15. In these plots the number of negative ions is n_ = 6. 
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And finally, in Figs. (12) and (13) we illustrate the effect of varying the number 
of impurity ions (concentrations) on the monomer distribution and the potential 
/(r). Here, again, we compare the 3D results with the ones obtained by the method 
of Section 6. Clearly, the agreement between the two approaches is good. 



30 



4*10"- 



3*10"- 



^ 2*10-5 



1*10 



■5 _ 



- 



^ y ^ 

,H-- / S . \ 

^.■■ *" \ * 

,+, / » . \ 

+ + + + +__+_ + > \»- 

.'^ ^ • I 

. _ -- • * • 1 

ID, n.=4900 4 Vi 

3D (60^), n.=4900 i \i 

■ ID, n.=400 "■■,'- 

3D (60^), n.=400 *\\ 

ID, n.=6 \".i 

3D (60^), n.=6 v,'. 

\'! 

-A I \ \ 

5 10 15 

r 

Figure 12: Monomer distribution for varying n_ 



20 



1.5 - 



■ ID , n.=6 

♦ 3D (60^), n.=6 
■--■ID, n.=400 

3D (60^), n.=400 
ID , n.=4900 

3D (60^ ), n.=4900 



' * *>. 



0.5 - 






10 



20 
r 



30 



40 



Figure 13: Electric Potential /(r) for varying n- 



31 



8 Discussion 

We have derived the mean-field equations for a coupled ionic-polymer system by 
doing a saddle-point analysis on a functional integral representing the partition 
function of the system, (11). This analysis shows that all mean- field level thermo- 
dynamical properties are obtained by extremizing an appropriate real-valued func- 
tional (14). Moreover, we have shown that the functional (14) possesses a single 
minimum, guaranteeing a unique solution of the coupled mean-field equations. We 
have also described two different numerical procedures for finding the mean-field 
solution, and have applied those to the problem of a charged polymer confined to a 
spherical cavity and immersed in an electrolyte solution. 

Although our calculations were intended as an illustration of the advantages of 
the approach, they have given us some interesting insight into the physics of con- 
fined polyelectrolytes. It has been suggested that materials consisting of spherical 
voids imbedded in a polymer gel can be used as so called "entropic trapping de- 
vices," in which macromolecules such as polymers and DNA could be trapped and 
separated |]22| , pSJ . Therefore, such materials may have important applications in 



specific biochemical trapping, or even as microreactors for applications in organic, 
bioengineering and combinatorial synthesis |2^ . 



It has been hypothesized |2^, |2^, |2^, |2^ that such trapping is a result of the 
higher conformational entropy "enjoyed" by the polymer in the large spherical void, 
as compared to the narrow channels connecting the voids within the gel. Our cal- 
culations reveal another important aspect of the trapping phenomenon — namely, 
from our results it becomes clear that electrostatic interactions also play a very 
important role in it. From Figs. (3), (5), (10) and (11) it is seen that polymers 
with high monomer charge p, or in dilute electrolyte solution, have a distribution 
function with a peak near the edge of the spherical cavity, while polymers with low 
monomer charge or in concentrated electrolyte solution (where the monomers are 
highly screened by the impurity ions, thus experiencing weaker repulsion from each 
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other), have a more flattened distribution, and are more hkely to be found near the 
center, rather than the edge of the cavity. In addition, highly charged monomers 
would possess higher energy in the voids, due to the stronger repulsion from each 
other as they are brought closer together in the folded polymer. Therefore, we 
expect that polymers with relatively low average monomer charge or in a concen- 
trated electrolyte solution would be easier to trap in spherical voids. Usually the 
monomer charge is approximately constant, and it is the impurity ion concentration 
that can be varied in the laboratory. We suggest that experiments be performed to 
investigate polymer trapping dependence on electrolyte concentration. 

We also expect that the lattice field theory approach to the statistical mechanics 
of charged polymers in electrolyte solution which has been developed in the present 
work will be useful for studying low symmetry systems involving complex liquids, 
such as colloidal suspensions in polyelectrolyte solutions. 



33 



References 

[1] R.D. Coalson and A. Duncan, J. Chem. Phys. 97, 5653 (1992). 
[2] A.M. Walsh and R.D. Coalson, J. Chem. Phys. 100, 1559 (1994). 
[3] N. Ben-Tal and R.D. Coalson, J. Chem. Phys. 101, 5148 (1994). 

[4] R.D. Coalson, A.M. Walsh, A. Duncan and N. Ben-Tal, J. Chem. Phys. 102, 

4584 (1995). 

[5] R.D. Coalson and A. Duncan, J. Phys. Chem. 100, 2612 (1996). 

[6] A.L. Kholodenko and A.L. Beyerlein, Phys. Rev. A 34, 3309 (1986). 

[7] A. Duncan and R. Mawhinney, Phys. Rev. D 43, 544 (1991). 

[8] R. Podgornik and B. Zeks, J. Chem. Soc. Faraday Trans. II 84, 611 (1988); R. 
Podgornik, J. Phys. A. 23, 275 (1990). 

[9] I. Borukhov, D. Andelman and H. Orland, "Random Polyelectrolytes and 
Polyampholytes in Solution" LANL preprint ( |cond-matt/9804304|) , 1998. 

[10] It was shown in [Q that an arbitrary dielectric profile e(r) can be incorporated 
into the LFT formalism in a straightforward way. For simplicity, we consider 
only the prototypical case of constant e here. 

[11] M. Doi and S.F. Edwards, The Theory of Polymer Dynamics (Oxford University 
Press, Oxford, 1986). 

[12] F.W. Wiegel, Introduction to Path-Integral Methods in Physics and Polymer 
Science (World Scientific, Singapore, 1986). 

[13] D.A. McQuarrie, Statistical Mechanics (Harper and Row, New York, 1976). 

[14] It is also possible to include a charge density term that accounts for one or many 
fixed macroions, such as the charged colloid particles studied in |l|]-[^. We will 



34 



not consider this complication in the present work, because the system upon 
which we focus, a charged polymer in an electrolyte solution, is interesting in 
its own right. 

[15] The condition of ground state dominance naturally requires that the polymer be 
confined by a potential (e.g., walls of a container), so that the energy spectrum 
of the Hamiltonian operator in ([13|) is discrete and M{Ei — Eq) ^ 1^ Ei being 
the energy eigenvalue corresponding to the first excited state. 

[16] G.H. Golub and CF. Van Loan, Matrix Computations (Johns Hopkins Univer- 
sity Press, Baltimore, 1996). 

[17] W.H. Press et al.. Numerical recipes in C: the art of scientific computing (Cam- 
bridge University Press, 1992). 

[18] Note that according to the notational conventions established in Section 2, Xc 
(cf. (|53D) is the negative of the standard electrical potential |]19|. Thus /(r) is 
large and positive near the negative charge source provided by the monomers 
of the polymer chain. 

[19] J. D. Jackson, Classical Electrodynamics, 2nd Ed. (John Wiley & Sons, New 
York, 1975). 

[20] S.A. Asher, J. Holtz, L. Liu and Z. Wu, J. Am. Chem. Soc, 116, 4997 (1994). 

[21] L. Liu, P. Li and S.A. Asher, J. Am. Chem. Soc, 119, 2729 (1997). 

[22] L. Liu, P. Li and S.A. Asher, Nature 397, 141 (1999). 

[23] L. Liu, P. Li and S.A. Asher, in press. 

[24] E.F. Casassa, Polymer Lett. 5, 773 (1967). 

[25] E.F. Casassa and Y. Tagami, Macromolecules 2, 14 (1967). 

[26] A. Baumgartner and M. Muthukumar, J. Chem. Phys. 87, 3082 (1987). 



35 



[27] M. Muthukumar and A. Baumgartner, Macroinolecules 22, 1937 (1989). 



36 



